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ABSTRACT 

A new method is presented for calculating the time evolution of spherically symmetric Type la 
Supernova in the post-explosion phase, enabling light curves and spectra to be simulated in a physically 
self-consistent way. The commonly exploited radiative equilibrium, that is in essence a gas energy 
balance condition, is unsuitable for this purpose for important physical and numerical reasons. Firstly, 
the RE depends on the heating and cooling rates of the gas by the radiation field, two quantities that 
almost completely cancel and are very difficult to calculate accurately. Secondly, the internal energy 
of the gas is only a tiny fraction of the total energy in the system (the vast majority of the energy 
resides in the radiation field), so that the vast majority of the energy is neglected in solving for the 
energy balance. The method presented in this paper, based on the radiation energy balance, addresses 
the bulk of the energy, does not depend on the heating/cooling rates, guarantees an accurate run 
of the bolomctric luminosity over time while bringing the gas temperatures into consistence with 
the radiation field. We have implemented the method in the stellar atmosphere code PHOENIX and 
applied it to the classical W7 model. The results illustrate the importance of each of the four physical 
contributions to the energy balance as a function of time. The simulated spectra and light curves for 
W7 show good resemblance to the observations, which demonstrates what can be done using PHOENIX 
with the REB method. 

Subject headings: stars: Supcrnovae: general - radiative transfer - methods: numerical 



1. INTRODUCTION 

Type la Supernovae (SNe la) can be divided physically 
into two phases: the explosion phase and the ballistic or 
free-expansion phase. 

In the first phase, a white dwarf explodes. During the 
explosion, the material in the white dwarf undergoes a 
series of thermonuclear reactions (Hoyle & Fowler (1960), 
Hillcbrandt & Nicmeycr (2000)). These reactions release 
a large amount of energy and produce heavy chemical 
elements. The energy that is released heats up the ma- 
terial and accelerates it to high velocities. The explosion 
process takes place in seconds. The resulting velocities 
are so large that the gravitational forces rapidly drop as 
the material expands, and the pressure gradients rapidly 
decrease as the material expands and cools. Thus the 
whole structure becomes ballistic, meaning that the ve- 
locity structure of the material no longer changes. 

In the second phase, another energy release mecha- 
nism becomes important. Instead of the expanding ejecta 
cooling down and becoming dim, the radioactive decay of 
unstable isotopes produced in the explosion phase heat 
it, making it glow (Truran et al. (1967), Colgate & Mc- 
Kee (1969)). The heating occurs via the gamma rays 
that are produced in the decay of these isotopes, which 
lose part of their energy as they pass through the ex- 
panding ejecta. The main contribution comes from 56 Ni, 
which decays through unstable 56 Co, to stable 56 Fe. 

It is the ballistic phase of SNe la that is observable, and 
it's luminosity is high enough to overpower that of the 
host galaxy. The extreme peak luminosities of SNe la, 
as well as their similarities in peak luminosity and spec- 
tra, and the correlation that exists between their peak 
luminosity and the rate they fade make SNe la very in- 



teresting objects observationally and theoretically. 

Successful theoretical accounts of the qualitative obser- 
vational features of SNe la date back to the early eight- 
ies (explosion phase: Nomoto et al. (1984) and ballistic 
phase: Branch et al. (1985)). Since then, the theory and 
methods for modeling both phases have evolved, but re- 
producing the detailed observational features of SNe la 
quantitatively is still a challenge. The impetus to in- 
crease the complexity and sophistication of the models 
will therefore continue for the foreseeable future. 

One recent development is the transition from ID to 
2D and 3D explosion models. In 3D explosion models 
using state-of-the-art mechanisms (like deflagration-to- 
detonation transition [see, e.g., Kasen et al. (2009)] and 
gravitationally confined detonation [see e.g., Jordan et al. 
(2008)]), significant deviations from spherical symmetry 
in composition and density are found. These angular 
inhomogencities require the radiative transport models 
for the ballistic phase also to be 3D, because it is not 
possible to determine the impact of a ID approximation 
(i.e., averaging over solid angle) without first doing the 
full 3D treatment. 

Although 3D radiation transport codes suitable for 
SN la models are available (Hauschildt & Baron (2010), 
Kasen et al. (2008)), the computational demands are cur- 
rently overwhelming. This limits the amount of detail 
that can be included in other areas; e.g., resolution in 
space, time and momentum, and the ability to treat 
the interactions between the radiation field and mat- 
ter micro-physically using statistical equilibrium (called 
non-local thermodynamic equilibrium (NLTE)), rather 
than using local thermodynamic equilibrium (LTE). The 
impact of such approximations on multi-D models can 
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be estimated from ID models, in which the physics other 
than the multi dimensionality can be treated in most de- 
tail. 

In subsequent papers, we will study the effects 
of NLTE including ionization that results from the 
non-Maxwcllian electron distribution. In this paper, 
we present a new method for solving the radiation- 
hydrodynamical energy equation (with hxed velocities), 
assuming spherical symmetry, that makes it possible 
to evolve SNe la from the end of the explosion phase 
throughout the ballistic phase, and therefore to simulate 
light curves and spectra, in a physically self-consistent 
way. 

2. OUTLINE OF THE PROBLEM 

SNe la in their ballistic phase are radiation- 
hydrodynamical systems. The hydrodynamical aspect 
becomes trivial if the velocities are fixed (which is a 
reasonable approximation 1 - hence the term "ballistic 
phase"). Given this approximation, the evolution of the 
density structure of the ejecta with time is then known a 
priori. Nevertheless, the fact that the material is moving 
is crucial in the treatment of the radiation. The dy- 
namical equations for a radiating fluid contain velocity- 
dependent terms of the order (v/c), unlike the situa- 
tion for non-radiating fluids, where the frame-dependent 
terms are only 0(v 2 /c 2 ) and are negligible even for the 
high velocities found in SN la ejecta, as pointed out by 
Castor (1972). 

Via the ideal gas law, and if we assume LTE, the caloric 
equation of state (i.e., the Saha-Boltzmann equations), 
the gas pressure and internal energy are related to the 
temperature, so that there is only one variable related to 
the properties of the gas left in the system 2 . There are 
no a priori constraints on the radiation field variables. 
The equations describing the system are the radiative 
transfer equation (RTE), the gas energy equation (GEE) 
and the radiation energy equation (REE). The mechan- 
ical energy equation, which is used to compute changes 
in the velocity structure, does not need to be considered 
in the ballistic phase since the velocities do not change. 

Although the general RTE explicitly includes tempo- 
ral variations of the radiation field, this term has been 
shown to often be unimportant in astrophysical flows 
(Castor (1972); Buchler (1979); Mihalas & Weibel Mi- 
halas (1984)) and in SNe la specifically by Baron et al. 
(1996). Kasen et al. (2006) found that SN la spectra ob- 
tained taking into account the full timc-dcpcndcnce do 
not differ much from time- independent calculations with 
the correct boundary conditions (inferred from the full 
calculations). Accordingly, the time-dependence of the 
RTE is weak (and neglected in this work), so that at 
every instant of time the radiation field is in large part 
determined by the temperature structure at that time. 

In this picture, the GEE and REE provide constraints 
on the temperature, involving the wavelength-integrated 

1 Pinto & Eastman (2000) and Woosley et al. (2007) studied 
post-explosion acceleration effects and find changes to the velocity 
structure of up to 10%. 

3 In NLTE the situation becomes more complicated, where for- 
mally tens of thousands of indirectly coupled gas variables (namely, 
the occupation numbers of all atomic levels) are present in the sys- 
tem. The method presented in this paper can naturally be applied 
to NLTE models, but we will postpone a detailed discussion of 
NLTE until the next paper (which is in preparation). 



angular moments of the radiation field. While the GEE 
and REE are time-dependent, the temperature at a sin- 
gle point in space-time is formally coupled, through the 
radiation field and in a highly non-linear way, to the tem- 
peratures at all other points in space and time before the 
current time. 

Therefore, a temperature-correction procedure is 
needed that converges to a solution of the temperature 
on the whole space-time grid that satisfies the GEE and 
REE conditions. 

One way to approach this problem was proposed in 
Jack et al. (2009) but it was not applied to a SN la 
problem. In their most realistic test scenario two of the 
four important energy terms (the flux gradient and the 
work done by the radiation field) together referred to as 
"structure term" were dropped (their energy equation 
22). In two follow-up papers (Jack et al. (2011, 2012)) a 
different approach, based on the GEE, is proposed and 
applied to a classical SN la problem. But the presented 
light curves strongly disagree with the observations in 
that the early magnitudes drop monotonically from the 
U-band all through the BVRIJHK-bands with 2.5 mag- 
nitudes (U-I Band in Figures 12-16 of Jack et al. (2011) 
and I-K bands in Figures 2-5 of Jack et al. (2012)), a fact 
that is obscured in the Jack et al. papers by shifting the 
observed light curve data to the model. Furthermore, 
in Jack et al. (2011) it is described that the method is 
too computationally intensive to calculate the whole light 
curve and that big forward jumps in time have to be 
taken in which the time-dependence (their Equation (1)) 
is not tracked. Consequently, the scheme is physically 
inconsistent and does not lead to accurate light curves 
and spectra. 

3. MODIFIED RADIATIVE EQUILIBRIUM: THE CAVEATS 

In the ballistic phase, there is only one unknown gas 
variable and that is the internal energy (or equivalcntly 
gas temperature or gas pressure), as explained above. 
Considering all mechanisms of energy transfer to and 
from the material (i.e., the first law of thermodynam- 
ics) leads to the gas energy equation. The GEE in the 
co-moving frame, accurate to 0(v/c), is given by (Miha- 
las & Weibel Mihalas (1984), Equation (96.7)) 



d d 1 

dt e+P dt-p. 



= 4nQ + pe 



(1) 



Here p is the density, e is the internal energy per unit 
mass, p is the gas pressure, e is the heating of the gas 
by gamma rays 4 per unit mass per unit time, and inQ 
is the net heating by radiation per unit volume per unit 
time, where 



(2) 



X\ and rj\ are the radiative monochromatic extinction 
coefficient and emissivity, J\ is the monochromatic mean 
intensity, and A is the wavelength. Equation (1) states 
that the rate of change of the material energy density 
plus the rate of work done by the material pressure equals 

5 The heating due to positron annihilation is included, and the 
gamma photons produced are transported like primary gamma 
photons. 
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Figure 1. The ratio of radiation energy density i? ra( j over gas 
energy density i? gaa is very high throughout a SN la structure at 
different times after the explosion (here the results for the W7 
model are shown). E gSjB oc pT does not gain as much from the 
high temperatures at low densities as the radiation energy density, 
which roughly scales like i? ra( j oc T 4 (see e.g. Mihalas & Weibel 
Mihalas (1984)). 

The low ratio values at low velocities reflect the absence of radioac- 
tive material around the center in the W7 model (mass fraction 
< 10~ 7 ), which at early times is too opaque for gamma rays to 
penetrate and thus remains cool. 

the net rate of energy input from the radiation field and 
thermonuclear sources. 

Evaluating the terms of Equation (1) in SN la models 
shows that after a few days from the explosion (depend- 
ing on the abundance of radioactive material throughout 
the structure) the terms on the left-hand side are negli- 
gibly small compared to the two terms on the right-hand 
side. The reason is that the internal energy density of 
the material is extremely small compared to the energy 
density of the radiation field, given the low densities and 
high temperatures in the ejecta, as shown in Figure 1. 
Thus to a good approximation (after the first few days) , 
the GEE (1) reduces to 



4ttQ = — pe 



(3) 



which is similar to but is not the strict radiative equilib- 
rium (RE) condition Q = 0. 

Equation (3) offers a tempting, simple way to deter- 
mine the temperature structure of the ejecta. The net 
heating term Q can be evaluated by doing the full radia- 
tion transport. Departures from the condition expressed 
by Equation (3) can be corrected for in a quasi-local 
fashion, since in zeroth order approximation, locally the 
net heating by the radiation field increases with decreas- 
ing local temperature. Corrections for non-local effects 
(e.g., the fact that the local net radiative heating effect 
decreases with decreasing global temperatures) can be 
accounted for iteratively. 

There are three major problems with this formal- 
ism. Firstly, the opacity changes dramatically over small 
wavelength scales (i.e., a line width), so that the wave- 
length sampling of the opacity becomes important in 
evaluating the integral in Equation (2). Also, the re- 
sult changes if lines are omitted. Secondly, because Q is 
a small difference between two large numbers [see Equa- 
tion (2)] - it is typically smaller by a factor of 1 3 — 10 7 
than the sum of both - the precise value of J near line 



centers becomes important, since there the weighting is 
the strongest. Therefore, this method is very sensitive 
to approximations in the radiative source function or in- 
accuracies in the solution of the RTE. Finally, the GEE 
is not the only energy balance condition. The vast ma- 
jority of the energy is stored in the radiation field, and 
this method does not take into account the local battery 
effect of the storage of energy in, and at later time the 
retrieval of energy from, the radiation field (see section 
6). 

Indeed, the gas can not absorb a significant amount 
of the energy deposited by radioactive decay, because 
through material-radiation interactions it is directly 
transferred to the radiation field, with little effect on 
the temperature. Let us therefore now turn our atten- 
tion from the material to the radiation, and consider the 
physical effects on the energy density of the radiation 
field. 

4. THE RADIATION ENERGY BALANCE (REB) METHOD 

4.1. The target flux 

The REE in the co-moving frame in spherical symme- 
try, accurate to 0(v/c), is given by (Mihalas & Weibel 
Mihalas (1984), Equation (96.8)) 



d J „d 1 .„_, T . v 

7T7- + ^ 7T7 (3A - ,/) — 

at p at p pr 



l__d 

r 2 dr 



(r 2 if), (4) 



where and J, H and K are the zeroth, first, and second 
angular moments of the radiation field, which physically 
represent the radiation energy density, the radial flux, 
and the radiation pressure (see e.g. Mihalas & Weibel 
Mihalas (1984)). Q can now (cquivalently) be inter- 
preted as the net cooling of the radiation field by the 
material. Equation (4) states that the rate of change of 
the radiation energy density plus the rate of work done 
by radiation pressure (the second and third terms on the 
left-hand side) equals the net rate of energy flowing into 
the radiation field from the material, minus the net rate 
of radiant energy flowing out of a fluid element by trans- 
port, all per unit mass. 

With the high expansion velocities typical of SN la, 
the radial coordinate r of a co-moving grid point quickly 
(within minutes or less) becomes proportional to its ve- 
locity 

r = vt , (5) 

where t is the time since explosion. The derivative of the 
inverse density can then be written as 



d 1 

dip 



3 1 

Tp 



(6) 



and the terms proportional to K on the left-hand side of 
Equation (4) cancel out. 

The next step is to express Q in the REE (4) in terms of 
e using the GEE (1). The terms on the left-hand side of 
the GEE are much smaller than the corresponding terms 
in the REE (the left-hand side) in the SN la problem (see 
figure 1) and can be dropped 6 . The net cooling by the 

6 Note that dropping the gas energy density terms is not at all 
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material now does not have to be evaluated (with all of 
the uncertainties in doing so we mentioned before) be- 
cause the material does not (significantly) absorb energy 
and therefore all of the energy flowing into the material 
flows out to the radiation field. 

Using Equations (3) and (6) and the spherical symme- 
try boundary condition H(r = 0) =0, we can integrate 
Equation (4) to obtain a target flux H tg (r) for all radial 
coordinates r on the grid, 



2 H ts (r) = 



pr l d J 
c dt p 



ct 



-J 



2 

r z p 

47T 



e dr 



(7) 



Through this target flux equation, the first angular mo- 
ment of the instantaneous local radiation field is related 
to the zeroth angular moment and its time derivative of 
the local and underlying shells. This is a highly indi- 
rect and non-local condition for the instantaneous local 
temperature, the more so since in a non-grey atmosphere 
the influence of the temperature on the radiation field in- 
volves a very large number of individual opacity sources. 

An example of the contribution of each of the terms of 
the REE in different regions of the atmosphere and for 
different evolution times after the explosion is given in 
Figure 2, where the W7 model Nomoto et al. (1984) was 
used for the explosion phase (see section 5 for details on 
the test setup). These results were obtained by solving 
the target flux equation using the method described in 
the next section. Figure 3 shows the evolution of the 
different terms in the REE, integrated over the whole 
structure, as a function of time. 

The target flux Equation (7) docs not take into ac- 
count the thermodynamic state at the time of the tran- 
sition from the explosion phase to the ballistic phase. In 
particular, the internal energy e of the gas at the time 
of transition determines the temperature of the structure 
before the local heating by gamma rays and/or the "ther- 
mal" radiation field takes over. In regions with abundant 
56 Ni, this happens only a few minutes after explosion, be- 
cause the cooling rate due to adiabatic expansion of the 
gas (which is proportional to e) is still very high. But in 
the outer shells that have no radioactive Nickel, it takes 
a few days before the ejecta are optically thin enough 
for these shells to be heated by gamma rays and thermal 
radiation originating in deeper shells. For such condi- 
tions (optically very thick outer shells at early times), 
the target flux Equation (7) cannot be expected to yield 
realistic results. However, this is not a problem for the 
radiation energy balance (REB) method because the con- 
tribution of these regions to the global energy balance is 
insignificant. 

The important point is that the energy released by 
early nuclear decays cannot escape from the ejecta right 
away; it is stored locally in the radiation field and re- 
leased later on. Thus the radiation field acts as a battery; 
this has an important effect on the temperature structure 
at later times, and therefore on the light curves and spec- 
tra. 

4.2. Translation to temperature corrections 

We do not presume to determine analytically the exact 
temperatures (or temperature corrections) that satisfy 

required for the REB method, but here supports focusing on the 
big energy pool, which is the radiation field. 
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Figure 2. This figure shows the contributions to the rates of en- 
ergy change of the four terms of the radiation energy equation (4) 
("Battery" (red) is the first and "Adiabatic" (green) the sum of 
the second and third terms on the left-hand side, while "Deposi- 
tion" (blue) is the first and "Flux" (black) the second term on the 
right-hand side) multiplied by r 2 , plotted against the velocity for 
the W7 model at 4, 8, 16, and 32 days after explosion (pe). These 
four terms correspond to the terms in the target flux equation (7) 
after taking the spatial derivative 8/ dr. The x-axis is limited to 
2 • 10 4 km/s, beyond that all terms are small. 
Note that the absolute rates decline rapidly with time. 
At early times, energy is stored in the radiation field, and the bat- 
tery term is negative. At day 4pe, the rate of storage of energy in 
the radiation field is about as high as the work term ( "Adiabatic" ) 
caused by the expansion of the co-moving frame (in the observer 
frame, work would be done on material moving through a radiation 
pressure gradient). These two effects together almost completely 
balance the heating ( "Deposition" ) of the radiation field by the gas 
from the energy deposited in the gas by gamma rays. As a result, 
the gradient of the radial flux ( "Flux" ) becomes very small. 
At later times, the battery term and the term become smaller due 
to cooling by expansion but continue to be important. The relative 
sizes of the terms are governed by the evolution of the temperature 
structure as a function of time. Methods that neglect these effects 
in the determination of the temperature structure miss important 
physics. 
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Figure 3. The upper panel shows the contributions of the four 
terms in the target flux equation (7) for the W7 model: the flux 
(black), the battery term (red), the adiabatic cooling of the radia- 
tion field (green), and the heating of the radiation field by gamma 
ray absorption (blue) respectively. See Figure 2 for the radial dis- 
tribution of these terms. Integration of the flux term (Figure 2) 
from the center outwards yields the total flux R 2 H(R) at the outer 
boundary R, which is the luminosity of the structure. After about 
55 days, the energy flowing into the ejecta from nuclear decays has 
gone into work and luminosity in equal amounts. 
Note that the bolometric luminosity (black curve) is a direct prod- 
uct of the REB method. 

The lower panel shows the integration over time of these same to- 
tal rates as a function of time (the quantities are assumed to be 
constant before t=l for simplicity). It shows that the radiation 
field battery term integrates to zero. This is an important verifi- 
cation of the radiation energy balance (REB) method. Physically, 
the battery term cannot release more energy than what has been 
stored previously, but this requirement is not directly accounted 
for, since the target flux equation is in differential form in t and 
not in integral form. Whereas instantaneous energy conservation 
(the requirement that the sum of the four rate terms be zero) is ex- 
plicitly solved for in each point in time, this result shows that the 
d(J/p)/dt computation is accurate and satisfies temporal energy 
conservation in this time-dependent term, individually. 



the target flux equation. Instead, we propose a method 
that translates departures of the actual flux from the 
target flux to approximate temperature corrections and 
iterate for convergence to the physically correct solution. 
The derivation partly follows Hauschildt ct al. (2003) 
and Lentz et al. (2003) where the Unsold-Lucy method 



(Unsold (1955), Lucy (1964)) is generalized to spherical 
geometry. 

We start from the observation mentioned before that 
the heating and cooling rates of the gas by the radiation 
are much bigger in size than their sum 



-»hcat 



-Q co ° l > Q 



heat 



= Q, 



(8) 



where <5 heat and Q c ° o1 represent the first and second term 
on the right-hand side of Equation (2) This condition 
does not contain enough physics to determine the dy- 
namics of the system (it contains even less physics than 
Equation (3)), but is adequate for the purpose of driving 
the model towards the physically justified condition of 
Equation (7). 

The extinction coefficient and cmissivity are assumed 
to consist of true absorption (k) and pure scattering (rr) 
contributions \\ = k\ + o~\, and thermal emission and 
scattering contributions rj\ = K\B\ + o~\J\, respectively, 
where B\ is the Planck function. Note that this macro- 
physical description is suitable for LTE but needs to be 
generalized for NLTE (see our paper on NLTE in prepa- 
ration). From Equations (2), (3), and (8) follows 



n\B\ dA 



K\J\d\ , 



(9) 



where the scattering contributions cancel. Let B be the 
wavelength-integrated Planck function. If one subtracts 
Equation (9) for the current temperature iteration from 
the next iteration (marked with primes) and assumes 
that the wavelength-averaged opacities, 

1 f°° 
= — / k\B x dA , 
B Jo 

1 f°° 
Kj = - I k\J\ dA , 
■J Jo 



Xh = jj J X\H\ dA , 

are insensitive to temperature changes (which is a good 
approximation, as observed by Lucy (1964)), Equation 
(9) can be written as 



AB=^AJ 

KB 



(10) 



where AB = B' — B and A J = J' — J. The temperature 
correction AT follows from AB using the temperature 
derivative of the wavelength integrated Planck law, 



AT= -J— AB= T -^ 
dB/dT 4 B 



"IT ■ (11) 



Following Unsold (1955), we derive AJ from the radi- 
ation momentum equation. In the co-moving frame and 
assuming spherical symmetry, this equation, accurate to 
0(v/c), is given by (Mihalas & Weibel Mihalas (1984), 
Equation (95.21)) 



d 
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XhH 



dv 
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H 



-U-l^iJ + K). (12) 
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Using the sphericity factor q introduced by Auer (1971), 



q = ~~f exp 



r 3/-l 
. fr' 



dr' 



where / = Kj J is the Eddington factor and r c is a ref- 
erence radius, the left-hand side of Equation (12) can be 
rewritten, and the right-hand side can be evaluated for 
Equation (5), yielding, 



J_d_ 

qr 2 dr 



(qfr'J) = -[XH + 



1 d 
H - --H 

cdt 



(13) 



Subtracting Equation (13) for the current temperature 
iteration from the next iteration and assuming that / is 
insensitive to the temperature one obtains, 



J_d_ 

qr 2 dr 



(qfr 2 AJ) 



XH 



ct 



AH 



(14) 



Here AH = H' - H = H tg - H. The time derivative in 
Equation (13), evaluated on a discrete time grid against 
the reference time t* , becomes 



At 



-H 



H te - H* 
t-t* 



H -H* 

t-t* 



AH 

t-t* 



(15) 



The numerator on the right-hand side is independent of 
the size of the time step in the denominator so that this 
term has to be dropped. This is the case because the 
reference value H* is the same for both temperature iter- 
ations. Integration of Equation (14) from the outermost 
model radius R to r yields 



AJ{r) 



1 



q(r)f(r) 



q(R)f(R)g(R) AH(R) 



- ^ f R l( r ') (xh(/) + |) r 2 AH(r') dr'j , (16) 

where g = J/H is the second Eddington factor, which is 
assumed to be insensitive to the temperature for the out- 
ermost model layer at radius R. From Equations (10), 
(11) and (16), we obtain the final expression for the tem- 
perature correction: 



AT = 



T K ) 



AB K B q(r)f(r) 



q(R)f(R)g(R) AH(R) 



(17) 



All of the quantities on the right-hand side of Equation 
(17) are available upon completion of each temperature 
correction iteration. Note that both H and H tE change 
with temperature, and iteration is required to converge 
H to H tg . 

The number of iterations required is found to range 
between 3 and 10, using the Active Damping method to 
accelerate convergence (paper in preparation), depend- 
ing on the assumed initial temperature structure and the 

7 Although this is typically a good approximation, it does not 
have to hold exactly, given the intentional approximate nature of 
the correction being calculated. 



relative contributions of the terms in the target flux equa- 
tion (the constant third term is the fastest to converge 
and the time-dependent first term is the slowest). 

4.3. Sequential solver for the time- dependence 

With the target flux equation (Equation (7)) and the 
temperature correction equation (Equation (17)), the 
formalism and the tools are available that are needed 
to solve, in principle, the radiation-hydrodynamical evo- 
lution of a SNc la during the ballistic phase. 

Two difficulties arise in the numerical evaluation of 
the time-dependence in the target flux equation (Equa- 
tion (7)). Firstly, the solution at location r and time t 
depends formally on the history of the whole structure, 
so that we face an initial condition problem. The initial 
condition for the time-dependence is 







(18) 



/=o 



which physically does not hold exactly, but is true to a 
very good approximation. 

Secondly, on a discrete time grid the time derivative 
can only be approximated, yet accuracy is important. 
Let x = J j p\ then the "both-sided" linear approximation 
is 

d Xi+i - Xi-i 



dt 



U 



(19) 



and the two "single-sided" linear approximations are 



■r, 



Xi-l 



dt' 
d 

dt fj-i_i — ti 



U — U-i 

Xi-\-\ Xi 



(20) 
(21) 



In the following, we refer to these three variants symbol- 
ically as >< , — > and <— , respectively. If the variation 
in J is non-linear, the both-sided approximation clearly 
is more accurate than the single-sided approximations. 
This is demonstrated in Figure 4. where both J/p and 
d(J/p)/dt are shown for the model described in sections 
5 and 6. If J progressively increases, — > under(over)- 
cstimates and <— over(under)-estimates the size of the 
positive slope, whereas if J progressively decreases, — > 
under (over)-estimates and 4— over(under)-estimates the 
size of the negative slope. 

We propose a sequential scheme to solve the time- 
dependence that makes use of the both-sided derivatives. 
The sequence of steps taken in the scheme are shown in 
the graph of Figure 5. There are two stages. 

In the first stage the proper startup conditions are de- 
termined, and in the second stage they are evolved to 
later times. The time derivative is set to zero initially 
and the temperature structure is computed individually 
for the time grid points. The problem of determining 
the proper startup conditions is equivalent to accurately 
determining the time derivative at early times after the 
explosion. As it turns out, the variation of J with time 
is highly non-linear at early times, so that relying on — > 
is systematically inaccurate. In order to reduce this sys- 
tematic effect, — > is replaced by the sequence — , 
where —X— acts on the starting and end grid points of the 
two — > operations, i = 1 and i = 3 respectively, yield- 
ing an improved value at point i = 2. In each of these 
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Figure 4. This figure shows the run of the radiation energy den- 
sity per mass over time for three different layers in the atmosphere 
on a normalized linear scale (lower panel). These curves have been 
computed with the model setup described in sections 5 and 6. 
Layer 95 is a rather deep, 56 Ni rich layer, layer 20 lays further 
outside, far from 56 Ni rich regions. Layer 50 does not contain 56 Ni 
but lies close to regions that do. In early times the structure is op- 
tically thick for gamma rays so that layers further out are not yet 
heated by the nuclear decay. Therefore, the rise in energy density 
sets in later in the outer layers 20 and 50 than in the deeper layer 
95. 

The slopes of these curves, determined using the 'both-sided' 
( >< ), and the two 'single-sided' (— > and <— ) approximations, are 
plotted in the upper panel on a logarithmic scale, where positive 
slopes are indicated with symbols and negative slopes without sym- 
bols. The single-sided approximations tend to over- or underesti- 
mate the slope in non-linear situations. 

three steps, the temperature structure is iterated to con- 
vergence. When the — > operation is used, the derivative 
is updated "on the fly" with the results of the updated 
temperature structure, Equation (20). With this result 
in hand, the derivative at i = 1 is updated with a — x — 
operation using the boundary condition (18). Whereas 
initially the temperature structure for i = 1 was derived 
assuming a time derivative of zero, equivalent to assum- 
ing a constant J/p over time, these steps need to be 
iterated until the derivative converges for i = 1. This 
typically takes about 4 cycles. 

In the second stage, the converged startup condition 
is evolved in time. Just as in the first stage, every — >■ 
operation is replaced by the sequence >< , so that 

every point is computed three times. The second and 
third time a point on the time grid is computed, the 
temperature converges quickly, because the corrections 
are relatively small (i.e., smaller than the changes from 
one to the next time step as in the first time). 
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Figure 5. This sequential scheme, consisting of two stages, de- 
scribes how the time-dependence is solved in the radiation energy 
balance method: (1) startup conditions (lower left part) and (2) 
the evolution of the startup conditions to later times (upper right 
part). The horizontal axis shows the time grid points i. Every dot 
in the graph represents a model step where the time derivative is 
"zero" , "single-sided" , or "both-sided" , which is indicated by the 
number of arrows pointing in, being zero, one, or two, respectively. 
The both-sided steps are indicated with thicker points. The star 
at i = is the fixed initial condition. The solid line shows the se- 
quence of model steps that need to be computed, while the dashed 
lines indicate the additional time-dependences. The dotted arrows 
indicate that the scheme can be trivially extended in that direc- 
tion. The dotted vertical line in the upper left corner means that 
the result provided by the lower left part is used for the first point 
of the upper right part of the scheme. 

Note that the use of <— has been avoided in the scheme. 
The reason is that when the <— operation is used, the 
derivative can not be updated on the fly with the results 
of the updated temperature structure, because the de- 
pendence of the target flux on the temperature is now 
inverted, which inevitably leads to poor convergence or 
even divergence of the temperature-correction method. 
Therefore, it would take much longer to achieve conver- 
gence using 4— operations. 

An alternative way to solve the time-dependence would 
be to solve the problem on the whole time grid simulta- 
neously. Again, the time derivative is set to zero initially, 
so that the temperature structure can be determined for 
every point in time independently Then at every point 
on the space-time grid, the time derivative is computed 
using — x— and the temperature structures are updated 
simultaneously This method is relatively computation- 
ally intensive, because in every temperature-correction 
iteration the radiation transport must be solved on the 
whole time (and space and momentum) grid, but it par- 
allelizes perfectly Therefore, the actual model runtime 
would not increase if the number of processors is scaled 
with the number of time grid points. Even though this 
method is less efficient, since it doesn't directly exploit 
the causal property of the evolution, it may be more 
accurate, because the use of single-sided steps can be 
completely avoided. Implementation of this alternative 
method has not been completed at this time, so we can- 
not give concrete results about its performance or bene- 
fits in accuracy. 

5. IMPLEMENTATION AND TEST MODEL ASSUMPTIONS 



We implemented the formalism and solution scheme 
described in the previous section in PHOENIX, a 
state-of-the-art, general purpose stellar atmosphere 
code (Hauschildt (1992), Hauschildt & Baron (1999), 
Hauschildt & Baron (2004), Hauschildt & Baron (2006)). 
Wide applicability of the code to astrophysical objects 
ranging from brown dwarfs to main sequence stars, AGB 
stars, novae and Supernovae, means the code has been 
well tested for a broad range of physical conditions. 
Specifically, a lot of work on SNe la has been done with 
PHOENIX, concentrating mainly on computing snapshots 
in time (e.g. Nugent et al. (1997), Lentz et al. (2001), 
Baron et al. (2006)) for explosion models like W7 Nomoto 
et al. (1984). PHOENIX solves the fully special relativistic 
radiation transport equation on a co-moving grid, includ- 
ing gamma-ray transport (see Lentz et al. (2001) and ref- 
erences therein), where a standard value for the effective 
absorptive gamma-ray opacity of k 7 = 0.06 Y e cm 2 g _1 
is assumed. Here Y e is the total electron number density 
divided by the baryon density. This method is less accu- 
rate than a full Monte Carlo treatment, but deviations 
are small as shown by Swartz et al. (1995). 

For the purpose of testing the REB method and quali- 
tatively investigating the effects of the different terms in 
the energy balance equation, we use the well known W7 
explosion model as presented in Branch et al. (1985). 
This model has become an established ID reference 
model that is known to reproduce qualitatively the ob- 
served spectral properties of SNe la. Furthermore, we 
assume LTE with a line absorption probability e = 0.3, 
based on Kascn et al. (2006), and defer the study of 
NLTE effects to a subsequent paper (in preparation). 
The ejecta are assumed to be in the ballistic phase, so 
velocities are constant in time. As mentioned before, we 
neglect the influence of the internal energy of the gas at 
the time of the transition from the explosion phase to the 
ballistic phase on the temperature in the early stages of 
the ballistic phase. 

We follow the evolution with time steps of 2 days, start- 
ing at t = 4 days. With time steps that are too big, the 
discrete approximations to the time derivatives (Equa- 
tions (19) to (21)) become inaccurate, while with time 
steps that are too small, the approximations become un- 
stable (because the denominator becomes small) . A time 
step of 2 days has been found to be a good compromise, 
but this choice may be optimized (e.g., allowed to vary 
with time) in the future to improve accuracy and perfor- 
mance. 

The number of points computed in the scheme (see 
Figure 5) is: 3 without time dependence, 17 in stage 1, 
and 69 in stage 2. The computation time in minutes per 
point in the scheme on 16 Opteron (2.2GHz) processors 
is: 30 for the first 3 points, 15 for the next 17 points, 
and 15 or 10 for the last 69 points (15 for the first time, 
10 for the second and third times a time grid point i 
is computed). The total time for calculating the light 
curves and spectra is thus 3 • 30 + 17 • 15 + 22 • 35 + 30) = 
1145 minutes, or 19 hours. Scaling beyond 16 processors 
is not very good 8 . Using 32 processors reduces the total 

8 For NLTE models the situation is different. These were found 
to scale well up to the number of layers used in the model, typically 
128. This is important, since NLTE models are computationally 
much more demanding. 



time by 15%. 

6. RESULTS 
6.1. Light curves 

Light curves and spectra for the W7 explosion model 
structure obtained using the REB method are shown in 
Figures 6 to 9 in comparison with observational data. 

In Figure 6 the model UBVRI light curves (we adopt 
the standard Bessel [1990] passband filters) are compared 
with MLCS light curve templates, as presented in Jha 
et al. (2007), and a close-up of BVR around maximum 
brightness showing the color properties is given in Figure 
7. The MLCS model assumes that the SN la light curve 
shape is correlated with the peak luminosity and can be 
described by the one-parameter function, 

Mjt(A) = M° X + P X A + Q X A 2 + 

5 - l0 ^( 65km/s°/Mpc ) ' ^ 

where is the template magnitude in photometric 

band X and H$ is the value of the Hubble constant. Mx 
depends only on the intrinsic luminosity-shape parame- 
ter A, and is a vector in the time domain. The vectors 
M. x , Px, and Qx are determined by a large sample of 
observations. Construction of the template is based on 
many observations , which provides a way to statistically 
minimize the uncertainties and the influence of the pe- 
culiarities of individual observations and objects. The 
templates describe, given the underlying approximation, 
the best observational knowledge about the typical shape 
and absolute magnitude of SN la light curves and their 
dependence on peak luminosity. 

The absolute observed magnitudes, provided by the 
templates, can directly be compared with the absolute 
magnitude obtained from the ballistic phase models. 
There arc no observational uncertainties like distance 
modulus or extinction in the templates. This means that, 
after aligning horizontally the time of observed maxi- 
mum brightness to the time of maximum brightness in 
the model light curves and choosing the maximum tem- 
plate magnitudes to correspond with the model using A, 
there are no free parameters that allow scaling of or off- 
sets in the magnitude, either in single bands or in all 
bands simultaneously. 

The obtained value for the luminosity-shape parameter 
is A = 0.08 ± 0.03 with H = 72 km/s/Mpc. This is 
a rather central value compared to the values for the 
observations used in the MLCS training, ranging from 
-0.40 to 1.38, where A < —0.15 is considered a "slow 
decliner" and A > 0.3 a "fast decliner" (see Jha et al. 
(2007)). 

The W7 light curves show good resemblance with the 
fitted MLCS light curves, meaning that the typical be- 
havior of objects with similar brightness is reproduced by 
the W7 model using the REB method. The most notable 
discrepancy is the late time behavior (after day 25pe) of 
the I-band; the synthetic flux in the I-band clearly is too 
large. 

6.2. Spectra 

In Figure 8, the model spectra are compared with 
observed spectra for SN 1994D (source: SUSPECT 
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Figure 6. Comparison of the synthetic light curves obtained for the W7 explosion model using the REB method with the MLCS template 
light curves (see also the close-up in Figure 7). The templates are based on a large sample of observations, in which each observation 
is individually corrected for reddening, extinction, etc. The one-parameter family of templates plus uncertainties represents the typical 
observed light curves and the sample variance in them for SNe la ranging from high to low peak luminosities. 

Direct comparison is made between synthetic and observational absolute magnitudes. The maximum B-band magnitudes are aligned in 
time and the MLCS parameter A = 0.08 ± 0.03 is fixed by matching the template and model maximum band magnitudes. Thus apart 
from t max and A, which are tightly constrained, there are no free parameters, no vertical shifts applied. 
The overall resemblance is good, especially given the absence of any free parameters in the fit. 

as 11 days before maximum B-band brightness, is little 
affected by dust extinction, has a high S/N ratio, and 
is often compared with W7 model spectra in the litera- 
ture. The observed fluxes have together been scaled with 
a single free parameter to account for the distance mod- 
ulus. Although not all of the features match well, the 
overall shapes of the observed spectra are reproduced 
well by the model spectra. Also, the luminosities of the 
spectra correspond well over the whole range in time, 
from 10 days before to 25 days after maximum B-band 
brightness. This indicates that the temporal evolution of 
the model luminosity accurately reproduces the observed 
evolution of the luminosity for this benchmark object. 

Single observations (and also single objects) can have 
peculiar features that may make comparison look better 
or worse than if other observations would have been cho- 
sen. Therefore, we also compare the synthetic spectra 
with the spectral templates from Hsiao et al. (2007) (see 
Figure 9). These templates are averages over many ob- 
servations, so that characteristic observational features 
dominate peculiarities. 

The synthetic spectra show a good resemblance to the 
templates both in terms of individual features and of the 
overall flux. Narrow features apparent in the synthetic 




Figure 7. This close-up of the B, V and R-band light curves 
of Figure 6 around maximum brightness shows that the colors of 
the model light curves (solid lines) reproduce the (observational) 
MLCS light curves (crosses) reasonably well. 

database). This was a nearby event, observed with very 
good spectroscopic temporal coverage starting as early 
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Figure 8. Comparison of the synthetic spectra obtained for the W7 explosion model using the REB method with observed spectra for 
SN 1994D. 

The date of the maximum B-band magnitude in the synthetic light curves [day 18 post-explosion (18pe)], is aligned with day after 
maximum B-band brightness (0pm). A single factor is used to scale all observed spectra simultaneously to match the synthetic fluxes. 
Apart from this factor there are no free parameters in this comparison. The good overall agreement of the fluxes shows that the evolution 
of the observed luminosity as a function of time is well reproduced using PHOENIX and the REB method. 
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Figure 9. Comparison of the synthetic spectra shown in Figure 8 with the spectral templates from Hsiao et al. (2007). The templates 
are averages over a large number of observations of SNe la, in order to let characteristic spectral features dominate over the peculiarities 
of single objects. 

A single factor is used to scale all of the templates simultaneously to match the synthetic fluxes. Apart from that there are no free 
parameters in this comparison. The good overall match of the fluxes shows that the evolution of the observed luminosity as a function of 
time is well reproduced using PHOENIX and the REB method. 
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spectra (e.g., between 4500 and 5000 in day 14pe) are 
missing or more "washed out" in the templates but arc 
present in the SN 1994D spectra in Figure 8. In the 
synthetic spectra the emission feature around 8500 is too 
strong on day 28pe and later. 

7. DISCUSSION AND CONCLUSIONS 

Up until recently, only snapshots in time could be done 
with Phoenix [see, e.g., Nugent et al. (1997)]. Such snap- 
shots are not physically self-consistent because they do 
not take into account the effects of the prior evolution, 
including the storage (and later release) of energy in the 
radiation field. Jack et al. (2011, 2012) propose a method 
to overcome these shortcomings but, while the method is 
based on the GEE and ignores the energy density of the 
radiation field, it suffers from serious physical and numer- 
ical problems as described in sections 2 and 3. The REB 
method proposed in this paper enables the PHOENIX ra- 
diation transport code - for the first time - to solve the 
time-evolution of Type la supernovae in an physically 
self-consistent and accurate way. 

The discrepancies in comparing the synthetic light 
curves and spectra calculated using PHOENIX and the 
REB method with observations are similar to those found 
with other codes (see, e.g., Kasen et al. (2006) and 
Kromer & Sim (2009)). In particular, the I-band is 
known to be very sensitive to the Ca II IR triplet (rest 
wavelength around 8580). A special treatment for the 
lines of this triplet can significantly improve the fit (see 
Kasen (2006)), the application of which is not limited to 
the W7 model (see, e.g., Woosley et al. (2007)). Note, 
however, that such tweaking is not based on physical ar- 
guments. Furthermore, NLTE effects may play a role 
as they are known to be important at later times. Yet, 
even though the W7 model is known to reproduce many 
observed properties reasonably well, it is a phenomeno- 
logical ID model lacking the detailed physics of many im- 
portant aspects of the of the explosion phase, and there 
is no reason to expect it to match the observations per- 
fectly. 

We have shown within the context of the W7 model the 
contributions made by the four physical processes that 
play a role in the evolution of the energy balance as a 
function of time. The contributions from the adiabatic 
cooling of the radiation field and the "battery effect" 
are large (much larger than the gradient of the radiative 
flux) and add at early times (<10pe). They then become 
smaller than the flux gradient and start to partially can- 
cel each other, but arc still significant (>10pe <60pe). 
Finally, they become so small that the flux transports 
away essentially all of the energy that is deposited lo- 
cally (>60pe). It is only after this time that the energy 
balance becomes time-independent because it's memory 
of the past fades away. Also it is only after about this 
time that the total amount of energy lost to radiative 
flux since day Ope becomes larger than the amount lost 
to adiabatic cooling through expansion of the radiation 
field. Neglecting these effects, as is done in snapshot cal- 
culations, will lead to incorrect temperature structures 
and therefore incorrect spectra. Furthermore, we have 
shown that the radiation-field battery term integrates to 
over time. This is an important verification of the ac- 
curacy of the REB method and our implementation of 
it, since the energy-balance equation at the heart of the 



method is in differential form in t, not in integral form. 

Generalizing the REB method to multi-dimensions is 
not straightforward. In multi-dimensions, there is a tan- 
gential flux in addition to the radial flux. Therefore, the 
target flux is no longer a scalar that can be obtained by 
radial integration, but a vector. 
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